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ABSTRACT 

These  lecture  notes  describe  a  variety  of  interesting  physical  phenomena  discovered  by  means  of  Direct 
Simulation  Monte  Carlo  (DSMC). 


1.0  DIRECT  SIMULATION  MONTE  CARLO 

Molecular  dynamics  (MD)  is  inefficient  for  simulating  dilute  gases  at  the  kinetic  scale.  The  main  reason  is 
that  the  relevant  time  scale  in  this  regime  is  the  mean  free  time  but  the  computational  time  step  in  MD  is 
limited  (by  numerical  stability)  to  the  collision  time.  Typically  there  are  many  orders  of  magnitude 
difference  between  these  time  scales  and  much  of  the  computational  effort  is  wasted  during  the  ballistic 
motion  between  collisions.  Molecular  dynamics,  on  the  other  hand,  is  efficient  for  liquids  since  molecules 
are  in  a  constant  state  of  interaction  with  their  surrounding  neighbour  molecules. 

Direct  Simulation  Monte  Carlo  (DSMC)  overcomes  the  inefficiency  of  MD  by  replacing  the  deterministic 
motion  with  a  stochastic  approximation  for  the  collision  process.  DSMC  is  able  to  advance  in  time  steps 
comparable  to  the  mean  free  time  between  collisions  yet  remain  accurate  at  the  level  of  the  Boltzmann 
equation.  Unlike  MD,  DSMC  is  always  numerically  stable  regardless  of  time  step. 

DSMC  was  developed  by  Graeme  Bird  in  the  late  1960’s  and  was  quickly  adopted  by  the  aerospace 
engineering  community  because  the  method  is  accurate  for  flows  with  high  Knudsen  number  (Kn),  the 
ratio  of  mean  free  path  to  system  length.  In  time  the  method  was  applied  to  an  expanded  number  of 
problems  in  physics,  chemistry,  and  engineering.  Examples  range  from  micron-scale  flows  (which  are  also 
high  Knudsen  number)  to  granular  gases  to  lunar  atmospheres.  The  algorithm  evolved  as  well,  with 
improvements  to  the  numerical  accuracy  and  efficiency  as  well  as  extensions  to  complex  chemistry  and 
even  to  dense  gases.  Nearly  50  years  after  its  introduction  DSMC  remains  the  dominant  numerical  method 
for  molecular  simulations  of  dilute  gases. 

The  details  of  the  DSMC  algorithm  are  presented  in  detail  by  Bird  [1]  and  in  tutorials  [2]  (as  well  as  in 
some  of  the  other  lectures  in  this  series)  but  for  the  sake  of  completeness  the  basic  scheme  is  outlined  in 
this  section.  DSMC  is  a  particle-based  scheme  so  a  typical  calculation  initializes  the  desired  geometry 
with  boundary  conditions  and  fills  the  computational  volume  with  random  particles.  At  each  time  step  all 
particles  move  ballistically  according  to  their  assigned  velocity;  in  DSMC  collisions  are  independent  of 
these  trajectories.  Any  particles  reaching  a  boundary  are  processed  according  to  the  imposed  conditions, 
such  as  randomly  assigning  a  new  velocity  to  a  particle  that  strikes  a  thermal  wall.  If  there  are  open 
boundaries  (e.g.,  wind-tunnel  configuration)  then  particles  are  generated  as  inflow  and  removed  if  they 
exit  through  the  boundary. 

The  core  of  DSMC  is  the  stochastic  (Monte  Carlo)  evaluation  of  the  collisions.  The  physical  domain  is 
partitioned  into  “collision  cells”  and  during  a  time  step  particles  are  randomly  selected  as  collision 
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partners  within  each  cell.  The  number  of  collision  to  occur  in  a  cell  is  determined  from  the  local  number 
density  and  temperature  (although  temperature  is  used  indirectly  through  average  relative  velocity).  A 
probability  that  a  pair  will  collide  is  assigned  based  on  their  relative  speed  and  collision  cross-section.  If 
the  collision  occurs  then  the  post-collision  velocities  of  the  colliding  particles  are  selected  at  random, 
preserving  conservation  of  momentum  and  energy.  If  the  particles  have  internal  degrees  of  freedom,  there 
are  further  details  as  to  the  re-assignment  of  energy;  chemical  reactions  add  further  steps  to  the 
computation. 

The  most  common  application  of  DSMC  is  in  aerospace  engineering  where  high  Knudsen  number  flows 
are  common  due  to  the  rarefied  gas  conditions.  Continuum  approaches  of  Computational  Fluid  Dynamics 
(CFD)  based  on  the  Navier-Stokes  equations  (and  its  variants)  are  not  accurate  for  these  so-called 
“transition”  regime  flows  because  the  stress  tensor  and  heat  flux  are  not  well  approximated  by  linear 
functions  of  gradients.  Another  common  application  of  DSMC  is  in  microscopic  flows  since  the  mean  free 
path  for  air  at  standard  conditions  is  roughly  0.05  microns.  An  example  of  this  type  of  flow  occurs  in  the 
lubrication  layer  between  the  read/write  head  and  the  spinning  platter  of  a  computer  disk  drive.  [3]  The 
head  is  lifted  from  the  platter  by  high  pressure  region  that  develops  between  them  as  the  air  is  dragged  by 
the  spinning  platter.  Because  the  sensitivity  of  the  magnetic  response  varies  exponentially  with  distance, 
the  head-platter  spacing  is  typically  less  than  20  nanometers  and  the  Knudsen  number  of  the  flow  is  order 
one. 

The  next  four  sections  describe  four  surprising  hydrodynamic  results  that  were  discovered  using  Direct 
Simulation  Monte  Carlo.  A  common  theme  is  that  the  author  has  worked  on  all  four  of  these  problems. 


2.0  ANOMALOUS  POISEUILLE  FLOW 

Poiseuille  flow  is  the  fluid  motion  in  a  confined  channel,  such  as  a  pipe,  driven  by  a  force  acting  parallel 
to  the  channel.  The  most  common  type  of  forcing  is  a  pressure  gradient,  with  the  flow  directed  from  high 
pressure  to  low.  However  we  will  instead  consider  Poiseuille  flow  driven  by  a  body  force,  such  as  a 
constant  gravitational  acceleration.  [4]  Furthermore,  we  consider  a  simple  rectangular  geometry  with 
thermal  walls  left  and  right  and  with  periodic  boundary  conditions  in  the  other  directions  (see  Figure  1). 


Figure  1:  Acceleration-driven  Poiseuille  flow;  flow  geometry  (left),  fluid  velocity  (right). 

The  acceleration  is  downward  so  in  this  configuration  the  fluid  exits  the  bottom  and  is  reintroduced  at  the 
top;  the  hydrodynamic  variables  (fluid  velocity,  pressure,  etc.)  only  vary  in  the  y-direction  (perpendicular 
to  the  walls)  at  the  steady  state.  This  state  is  reached  when  the  momentum  added  by  the  body  force  is 
balanced  by  the  momentum  removed  by  the  viscous  drag  of  the  walls.  Similarly,  at  the  steady  state  the 
viscous  heat  produced  by  the  shearing  of  the  fluid  is  balanced  by  the  cooling  of  the  walls. 
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From  Figure  1  we  see  that  the  profile  of  fluid  velocity  measured  in  DSMC  is  in  good  qualitative 
agreement  with  the  quadratic  profile  predicted  by  Navier-Stokes  (once  the  slip  velocity  at  the  boundary  is 
accounted  for)  even  though  the  distance  between  the  walls  is  only  10  mean  free  paths  (Kn  =  0.1).  However 
the  temperature  profile  (see  Figure  2)  is  qualitatively  different  in  that  there  is  a  dip  in  the  centre.  Since 
heat  is  generate  within  the  gas  due  to  shear  and  can  only  be  removed  at  the  walls  this  means  that  near  the 
centre  of  the  channel  heat  is  flowing  against  the  temperature  gradient  (i.e.,  flowing  from  cold  to  hot). 


Figure  2:  Acceleration-driven  Poiseuille  flow;  temperature  (left),  pressure  (right). 

The  pressure  profile  is  also  anomalous  in  that  there  is  a  pressure  gradient  towards  the  walls  yet  at  the 
steady  state  there  cannot  be  any  fluid  velocity  in  this  direction.  The  Navier-Stokes  equations  have  no  term 
in  the  momentum  equation  to  balance  this  pressure  gradient  however  at  the  Burnett  level  we  find 
agreement  with  the  DSMC  measurements.  [5]  For  the  temperature  profile  one  has  to  go  to  even  higher 
order,  namely  super-Burnett  theory,  to  obtain  the  correct  profile  [6];  the  central  dip  in  temperature  is  also 
predicted  from  BGK  theory  [7].  Finally,  pressure-driven  driven  Poiseuille  flow  is  more  complex  because 
the  flow  is  not  one-dimensional  however  similar  effects  are  observed  [8];  the  profiles  of  temperature  and 
pressure  disagree  with  Navier-Stokes  predictions  but  are  in  good  agreement  with  super-Burnett  theory  [6]. 


3.0  ANOMALOUS  COUETTE  FLOW 

Couette  flow  is  the  flow  in  a  confined  channel  that  is  driven  by  a  shearing  of  the  fluid  created  by  the 
motion  of  the  channel  walls.  In  laboratory  experiments  this  is  most  easily  achieved  between  concentric 
cylinders  with  the  inner  cylinder  rotating  at  a  constant  angular  speed  (see  Figure  3).  There  are  a  number  of 
interesting  hydrodynamic  instabilities  that  occur  in  cylindrical  Couette  flow  however  we  will  restrict  our 
attention  to  low  Reynolds  numbers,  well  below  any  instability.  As  in  the  previous  section  we  consider 
high  Knudsen  numbers  (Kn  >0.1)  and  investigate  the  resulting  velocity  profile  of  the  gas. 

The  velocity  of  a  gas  moving  past  a  solid  surface  has  a  velocity  at  that  surface  that  differs  from  the 
surface’s  velocity,  even  if  the  solid  fully  thermalizes  molecules  reflecting  from  it.  This  effect  was 
predicted  by  Maxwell  and  confirmed  by  Knudsen;  the  physical  origin  is  the  difference  between  the 
impinging  and  reflected  velocity  distributions  of  the  gas  molecules  due  to  gradients.  We  may  define  a  slip 
length  (see  Figure  3),  which  is  the  distance  within  the  wall  at  which  the  velocity  profile  of  the  gas 
extrapolates  to  match  the  wall’s  velocity.  For  a  fully  thermalizing  surface  the  slip  length  is  approximately 
one  mean  free  path  so  it  is  only  significant  for  high  Kn  flows.  Slip  increases  if  some  particle  reflect 
specularly  and  we  define  the  accommodation  coefficient,  a,  for  a  surface  as  the  fraction  of  thermalized 
(non-specular)  reflections  from  that  surface.  [9] 
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Figure  3:  Cylindrical  Couette  flow  geometry  (left);  Slip  length  for  a  stationary  thermal  wall  (right). 


DSMC  simulations  of  cylindrical  Couette  flow  confirm  that  when  the  accommodation  coefficient  is  one 
(fully  thermalizing)  the  fluid  velocity  is  greatest  near  the  inner  (moving)  cylinder  and  decreases  with 
increasing  distance  from  the  center  (see  Figure  4).  On  the  other  hand  when  a  <  0.1  the  angular  speed 
increases  with  increasing  radial  distance;  when  the  walls  are  nearly  completely  specular  (a  =  0.01)  the  gas 
rotates  like  a  solid  body  (v  =  cor).  This  result  is  not  anomalous;  it  was  predicted  by  Maxwell  over  a  century 
ago  as  a  consequence  of  the  centrifugal  force  acting  as  a  body  force  on  the  fluid. 


Diffusive  Waite 


Specular  Walls 


Figure  4:  Couette  flow  regimes  (left);  fluid  velocity  versus  radius  (right). 


The  unexpected  result  discovered  by  DSMC  is  that  for  a  certain  range  of  parameters  the  fluid  velocity  has 
a  minimum  speed  within  the  fluid  (see  Figure  5).  [10]  This  effect  occurs  when  the  walls  are  80%-90% 
specular  and  has  been  rigorously  confirmed  using  BGK  theory.  [1 1]  At  present  there  is  no  simple  physical 
interpretation  for  this  phenomenon,  which  involves  the  interplay  of  the  curved  geometry  and  the 
asymmetric  velocity  distribution  in  the  gas.  Yet  it  highlights  the  unusual  and  unexpected  behaviour  that 
can  arise  in  transition  flows. 
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Figure  5:  Anomalous  Couette  flow  (left);  fluid  velocity  versus  radius  (right). 


3.0  ANOMALOUS  TEMPERATURE  GRADIENT  FLOW 

In  DSMC,  as  in  all  molecular  simulations,  the  computation  calculates  the  positions  and  velocities  of 
particles  as  functions  of  time.  But  the  variables  of  interest  are  hydrodynamic  quantities,  such  as  pressure 
and  temperature.  For  some  quantities  the  transformation  is  unambiguous,  for  example  the  mass  density  is 
the  total  mass  of  particles  in  a  local  region  divided  by  the  region’s  volume.  Other  quantities,  such  as 
temperature,  entropy,  Gibbs  free  energy,  etc.,  are  less  obvious.  Even  fluid  velocity  has  two  interpretations. 

The  intuitive  measure  of  fluid  velocity  is  the  average  velocity  of  particles  in  a  local  region  (i.e.,  a  cell), 
which  we  may  write  as, 


v  = 


v. 


where  the  sum  is  over  particles  within  the  cell.  The  same  result  is  obtained  from  the  centre  of  mass 
velocity  in  the  cell, 


J 

u  =  — 
M 


N 

zeC _ 

mN 


These  two  equivalent  expressions  give  the  instantaneous  fluid  velocity  in  the  cell;  the  mean  value  over  S 
samples  (either  time-averages  or  ensemble-averaged)  may  be  written  as, 


~S 


Yet  this  is  not  the  only  way  to  define  the  mean  fluid  velocity.  In  fact  it  is  not  the  definition  normally  used 
in  DSMC.  Instead,  the  mean  fluid  velocity  is  calculated  as  a  cumulative  average, 
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These  expressions  for  the  mean  fluid  velocity  are  not  equivalent  because  the  former  is  the  average  of  a 
ratio  (<momentum/mass>)  and  the  latter  is  the  ratio  of  the  averages  (<momentum>/<mass>). 

Consider  a  closed  box  filled  with  a  dilute  gas  and  subjected  to  a  constant  temperature  gradient  imposed  by 
thermal  walls.  Using  the  cumulative  mean  to  define  fluid  velocity,  DSMC  calculations  yield  the  expected 
result:  no  flow  at  the  steady  state.  However,  the  mean  instantaneous  fluid  velocity  gives  an  anomalous 
result:  The  fluid  has  a  velocity  that  rises  quadratically  to  a  maximum  in  the  centre  of  the  system  and  the 
magnitude  of  this  flow  varies  linearly  with  the  temperature  gradient  (see  Figure  6).  [12] 


Figure  6:  Anomalous  flow  in  a  temperature  gradient.  Mean  instantaneous 
velocity  (left);  Cumulative  mean  fluid  velocity  (right). 


The  physical  origin  of  this  effect  is  that  out  of  equilibrium,  as  in  a  temperature  gradient,  the  fluctuations  of 
mass  density  and  of  momentum  are  correlated.  At  equilibrium,  density,  fluid  velocity,  and  temperature  are 
conjugate  hydrodynamic  variables;  when  measured  simultaneously  they  are  uncorrelated.  However,  out  of 
equilibrium  there  are  asymmetries  that  produce  correlations,  which  can  be  computed  using  Landau- 
Lifshitz  fluctuating  hydrodynamics  and  that  are  confirmed  by  DSMC  simulations  [13].  For  example,  in  a 
temperature  gradient  the  density  fluctuates  above  average  when  the  fluid  velocity  fluctuates  in  the 
direction  against  the  gradient  with;  specifically,  density  and  fluid  velocity  fluctuations  are  correlated  as, 


(5p(x)5u(x))  oc  -x(L-x)VT 

The  two  definitions  of  fluid  velocity  may  be  related  as, 


1  + 


5N7 


2\\ 
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{ SJSN )  


ML  - 
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~w 


DSMC  simulation  results  verify  that  the  anomalous  fluid  velocity  is  entirely  explained  by  this  non¬ 
equilibrium  correlation  of  fluctuations. 

Finally,  it  should  be  noted  that  a  similar  bias  occurs  when  using  any  instantaneous  hydrodynamic  quantity, 
such  as  temperature  or  pressure.  [14]  A  similar  bias  is  also  created  if  boundary  conditions,  such  as  inflow 
from  a  reservoir,  do  not  correctly  include  thermodynamic  fluctuations.  [15]  Because  the  error  goes  as  1/A, 
where  N  is  the  number  of  particles  in  the  sampling  cell,  for  engineering  applications  it  may  be  acceptable 
compared  with  other  approximations.  However,  it  is  an  effect  that  all  DSMC  users  should  be  wary  of;  we 
first  thought  it  was  a  bug  in  the  code! 


4.0  ANOMALOUS  DIFFUSON  FLOW 

As  discussed  in  the  previous  section,  fluctuations  are  enhanced  when  a  system  is  out  of  equilibrium,  such 
as  when  a  gradient  is  imposed  by  boundary  conditions.  Figure  7  illustrates  this  effect,  showing  snap-shots 
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for  an  equilibrium  system  with  a  concentration  gradient  induced  by  gravity  and  a  non-equilibrium  system 
with  a  comparable  concentration  gradient  imposed  by  boundary  conditions.  Clearly  the  fluctuations  are 
significantly  greater  in  the  non-equilibrium  steady-state  scenario.  The  effect  of  fluctuations  is  even  more 
noticeable  in  non-steady  problems,  such  as  the  mixing  of  initially  segregated  gases  (see  Figure  8).  These 
giant  fluctuations  in  diffusive  mixing  have  been  observed  in  laboratory  experiments,  creating  variations 
visible  to  the  naked  eye.  [16] 


Figure  7:  Equilibrium  concentration  gradient  induced  by  gravity  (above); 
Steady-state  concentration  gradient  imposed  by  boundary  conditions  (below) 


Figure  8:  Snapshots  of  the  concentration  in  the  diffusive  mixing  of  two  gases  (red  and  blue)  at 
t  =  1,  4, 10  (top,  middle,  bottom),  starting  from  a  flat  interface  (phase-separated  system)  at  t  =  0. 


An  unexpected  consequence  of  these  non-equilibrium  fluctuations  is  that  they  enhance  the  effective  rate  of 
diffusion.  The  simplest  case  is  “red-blue”  diffusion  in  which  the  two  species  are  physically  identical.  In 
this  case  the  fluctuations  of  fluid  velocity  are  the  same  as  in  equilibrium.  Formulated  in  terms  of 
fluctuating  hydrodynamics  the  concentration  equation  is  coupled  to  the  velocity;  in  the  isothermal, 
incompressible  approximation  the  correlation  of  concentration-velocity  fluctuations  is,  in  Fourier  space, 


oc  — 


k2±Vc 

~ 


The  total  mass  flux  for  concentration  may  be  written  in  the  form,  [17] 


(y)  ~  (D0  +  AD)Vc  where  AD  = - ^—^(Scdu^dk 

8 n  k 

We  see  that  there  are  two  contributions  to  this  flux:  the  “bare”  diffusion  coefficient  D0  and  the 
contribution  due  to  the  correlation  of  concentration  and  velocity  fluctuations.  Because  the  latter  is  also 
linear  in  the  concentration  gradient  it  may  be  written  as  an  enhancement  to  the  diffusion  coefficient,  AD. 
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For  a  slab  geometry  (Lz  «  Lx  «  Ly )  we  have  that  A D  oc  In  Lx.  That  is,  as  we  make  the  system  wider  in  a 
direction  perpendicular  to  the  gradient  the  diffusion  is  enhanced;  A D  increases  because  the  integral  over 
wavenumber  extends  to  smaller  values  of  k  as  the  system  width  increases.  This  enhanced  diffusion  is 
anomalous  since,  in  the  deterministic  hydrodynamics,  the  diffusion  coefficient  is  independent  of  the 
system  geometry. 

This  effect  is  seen  in  DSMC  simulations  where  we  may  separate  the  contributions  to  the  concentration 
flux.  Specifically,  in  DSMC  we  may  independently  measure  the  averages  and  correlations  of  species 
density  and  velocity,  as  well  as  total  mass  flux  of  each  species.  As  seen  in  Figure  9,  the  effective  diffusion 
coefficient  increases  in  qualitative  agreement  with  theory  up  to  system  widths  approaching  the  system 
height  (Lx  ~  Ly).  [17] 


Figure  9:  Simulation  geometry  (left).  DSMC  measurements  of  bare  and  total 
effective  diffusion  coefficients  as  a  function  of  system  width  (right) 


5.0  CONCLUDING  REMARKS 

Direct  Simulation  Monte  Carlo  is  an  excellent  numerical  tool  for  both  applications  and  for  basic  research. 
Given  DSMC’s  computational  efficiency  even  modest  resources  are  adequate;  except  for  the  results  in 
Figure  9  all  of  the  results  presented  in  these  notes  may  be  reproduced  in  a  day  or  two  using  a  laptop. 
Furthermore,  DSMC  is  relatively  simple  to  program  and  modify;  public-domain  source  codes  are  available 
at  a  number  of  sites  (including  the  author’s).  Finally,  DSMC  is  an  excellent  educational  tool  for  teaching 
kinetic  theory  since  the  algorithm  is  as  simple  as  molecular  dynamics  yet  is  efficient  enough  to  yield  basic 
results,  such  as  the  measurement  of  viscosity  for  a  hard  sphere  gas,  in  a  matter  of  minutes. 
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